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ABSTRACT 

NBody realizations of axisymmetric collisional galaxy cores (e.g. M32, M33, NGC205, 
Milky Way) with embedded growing black holes are presented. Stars which approach 
the disruption sphere are disrupted and accreted to the black hole. We measure the 
zone of influence of the black hole and disruption rates in relaxation time scales. We 
show that secular gravitational instabilities dominate the initial core dynamics, while 
the black hole is small and growing due to consumption of stars. Later, the black hole 
potential dominates the core, and loss cone theory can be applied. Our simulations 
show that central rotation in galaxies can not be neglected for relaxed systems, and 
compare and discuss our results with the standard theory of spherically symmetric 
systems. 
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1 INTRODUCTION 

Galaxy cores are the hosts of supermassive black holes 
(SMBHs), the engines of quasars and of active galactic 
nuclei. There is increasing evidence that SMBHs play an 
important role in the formation and global evolution of 
galaxies. They are co mmonly observe d at the centers of 
many nearby galaxies (|Shankarl [20091 ') . and the existence 
of quasars at least up to redshifts z = 6 (|Degraf et al.l 
l20ld : IWillot et alj|20ich implies that many of these SMBHs 
reached nearly their current masses at very early times. 
Evolution of galactic nuclei during and after the era of 
peak quasar activity therefore took place with the SMBHs 
already in place. The energy released by SMBHs during 
and after the quasar epoch must have had a major im- 
pact on how gas cooled to fo rm galaxies and galaxy clus- 
ters (|Scannapieco et all 120051 ). However, the detailed his- 
tory of SMBH growth is still being debated. Some work 
has focused on the possibility that the seeds of SMBHs 
were black holes of much smaller mass — either remnants 
of the first generation of stars, so-called "Population III 
black holes" (|Madau fc Reedl200lh . or the (still speculative) 
"intermediate-mass black holes" (IMBHs), remnants of mas- 
sive stars that form in dense clusters via physical collisions 
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betw een stars (|Portegies Zwart et al.l |2004| ; IMapelli et al.l 
l2010h . 

Observations with the Hubble Space Telescope have elu- 
cidated t he run of stellar dens i ty near the centers of nearby 
galaxies (IFerrarese et al. I l2006l ; ICdte et al.ll2007l ; iGlass et all 
2011). Nevertheless, in the majority of galaxies massive 
enough to contain SMBHs, the central relaxation time is 
much greater than the age of the universe, due both to 
the (relatively) low stellar densities and also to the pres- 
ence of a SM BH, w hich increases « rms (|Faber et al 1 1l997l : 
IFerrarese et alj|2006l ) . This long relaxation times imply that 
nuclear structure will still reflect the details of the nuclear 
formation process. Beyond the Local Group, essentially all 
of the galaxies for which the SMBH's influence radius is spa- 
tially resolved have 'collisionless' (non-relaxed) nuclei with 
low nuclear densities. The nuclei of these 'core' galaxies may 
have been much denser before the cores were created by 
probably binary SMBHs. 

Only the smallest galaxies known to harbor SMBHs 
have nuclear relaxation times of ~ 10 Gyr or shorter. In 
such a nucleus the stellar distribution will have had time to 
evolve to a collisionally relaxed system. Galactic spheroids 
fainter than Mv — -18 show this property. The Milky Way 
nucleus is 'collisional'. It has a half- mass relaxation time 
t r h ~ 5 x 10 10 yr but the steep density profile implies 
t r h ~ 6 x 10 9 yr at 0.2 r^m (0.6 pc), where rh m is the 



2 J. Fiestas, O. Porth, P. Berczik, R. Spurzem 



half- mass radius, and ~ 3. 5 x 10 9 yr at .1 r^m (0.3 pc), 
assuming Solar-mass stars (jMerritt llioO^ ). Collisional nu- 
clei are present in th ree other Local Group galaxies (M32 , 
M33 and NGC205) (|Lauer et all 1 19981 ; IValluri et all I2005T I 
although M32 is the o nly one of these to exhibit dynamical 
evidence for a SMBH (jValluri et al.ll2005l l. 

Since the 1980s, the dominant model for the formation 
of elliptical galaxies and bulges - the stellar syste ms that 
conta in SMBHs - has been the merger model (|Toomrel 
1 19771 ). An almost certain consequence of a merger is 
the in-fall of the progenitor galaxies' SMBHs into the 
nucleus of the merge d system, resulting in the formation 
of a binary SMBH ijBeeel man et aD Il980f). as observa- 
tions of uncoalesced dual SMBHs show (|Rodriguez et al.l 
120061 ; IValtonen et~aH I2008T ). There is no generally agreed 
idea on whether and how fast a binary SMBH coalesces 
after a galaxy merger. Earlier work (mostly numerical 
simulations) in spherically symmetric nuclei discussed a 
stalling problem which would hang up the binary SMBH 
at some sub -pc separation (stal li ng or last parsec prob- 
lem, c f. e. g. iGould fc Rixl (|2000|); IMilosavlievic fc MerrittJ 
(l2003h; [Hemsendorf. Sigurdsson. fc Spurzeml d2002f); 
iMakino fc Funatol (|2004 ); iBerczik. Merritt. fc Spurzeml 



( 20051 )). But it turns out that any degree of more re- 



alism helps clearing the last parsec prob lem, such as 
axisymmetry of the galaxy merge r remnant (IBerczik et al.l 



20061; iPerets fc Alexander! |200S| ; iBerentzen et all l2009t 



Berczik et all 120111 ) or the presence of gaseous mate- 



rial in the nu c leus |Dotti et all 120091 ; iMaver et ail |2010| ; 
ICallegari et al. | 120091 . |201Q| ~ A self-consistent study of 
the combined effect of high-accuracy high-resolution 
stellar dynamics with binary black holes, together with 
evolution of a central galactic nuclear disk, and its 
interaction with stars and black holes is still lack- 
ing to our knowledge. Steps towards this goal could 
be in our view the p i oneeri ng semi-analytical study 
of IVilkoviskii fc Czernvl (|2002|) or the detailed gas- 



dynamical of 



and ste llar 

(|2009l ); IJohansson. Naab. 
latter lack the proper resolution 



Naabl 

( 20091 ), but the 
of previously cited 
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purely stellar dynamical work to follow the sub-parsec 
evolution of the binary SMBH well. One of the most 
detailed descriptions of physical processes between 
stars and g a s in galacti c nuclei has been prese n ted b y 
Ciotti et all (|2009|, |2010|); IShin. Ostriker fc Ciottll (|2010l) ; 
Shin. Qstriker.fc Ciottil (|2010l ). Their central resolution 



is mu c h better than that of Johansson. Burkert. fc Naabl 
[|2009h ; |jonansson. Naab.fc Burkertl (|2009h . but their mod- 
els are spherically symmetric and unable to resolve properly 
the collisional evolution of binary SMBH in dense gas-star 
sys tems in ga la ctic n uclei. Another preliminary approach 
of lJust et al.1 <|201lD follows the star-disk interactions 
with direct high-accuracy stellar d ynam ical models (as a 
follow up to IVilkoviskii fc Czernvl |2002l ). but still uses a 
stationary disk model). 

In this study we investigate following mechanisms im- 
portant in determining the structure and evolution of colli- 
sional galactic nuclei embedding growing SMBHs 



• Destruction of stars by the SMBH. A SMBH defines 
a "loss cone" of orbits that pass within its event horizon 
or tidal disruption sphere, r ^ rv Indeed the existence 



of a capture sphere is crucial for solutions like Bahcall & 
Wolf's, since it precludes the formation of an isothermal 
(/ ~ e E ) distribution of velocities which w ould necessaril y 
have a very high density near the SMBH |Peebleslll972h . 
Continued loss of stars to the SMBH also i mplies that no 
precisely steady-state equilibrium can exist l|Shapirolll977l ; 
iBaumgardt et al.ll2006l ); the nucleus will expand, on a relax- 
ation time scale, due to the effective heat input as stars are 
destroyed (|Shapirolll977l ). 

Continued supply of stars to the SMBH requires some 
mechanism for loss cone re-population. The most widely dis- 
cussed mechanism is gravitational encounters, which drive 
a diffusion in energy (E) and ang ular momentum (J) . 
The latter dominates the loss rate (|Frank fc Reesl Il976l ; 
iLightman fc Shapirolll977l ). Roughly speaking, many of the 
stars dominated by the BH potential (i.e. within its in- 
fluence sphere) will be deflected into r t in one relaxation 
time, i.e. the loss rate is roughly (M./m*)(t r h) _1 . In a colli- 
sional nucleus with M. « 10 6 M Q , this is ~ 10 6 /(10 10 yr) w 
10~"V -1 - 

Stellar disruption has direct observational consequences. 
Tidally disrupted stars produce X- and UV radiation with 
luminosities of ~ 10 ergs -1 , potentially outsh ining their 
host galaxies for a period of days or weeks (|Kobavashil 
|2004 Khokhlov fc Merle] 1 19961 ). A handful of x-ray flar- 
ing eve nts have been observe d that have the expected sig- 
nature l|Komossa et al.l [2004') , and the number of detec- 
tions is cru dely consistent with th eoretical estimates of the 
event rate (|Wang fc Merr itt 2004). Tidal flaring events may 
dominate t he x-ray luminosity function of AGN at Lx < 
10 44 ergs _1 (|Reeslll988l ; IMilosavlievic et al.ll2006l ). Compact 
objects (neutron stars or stellar-mass black holes) can re- 
main intact at much smaller distances from the SMBH; 
these objects would emit gravitational waves at poten- 
tially observable amplitudes before spiralling in and may 
dominate the event rate for low-frequency gravitational 
wave interferometers like LISA (|Hopman fc Alexanderll2006l ; 
lEilon. Kupi fc Alexande3l2009h . 

That loss cones can be much more quickly refilled in 
even only slightly axisymmetric or triaxial nuclei had 
been realized much earlier in the c ontext of tidal accre- 
tion of stars onto single black holes (iNorman fc SiTklll983l ; 
lYu fc Tremainell2002l ; lMerritt fc Poonll2004l ) . It is quite nat- 



ural after a galaxy merger that the merger remnant is 
not spherically symmetric and is not completely void of 
gas, so one would expect frequent mergers of supermassive 
black holes as well. It seems this is consistent with the cos- 
mological evolution of g alaxy and black hole populations 
(|Hirschmann et al]|2010l ). and also the LISA gravitational 
wave community expe cts frequent coalescences of binary 
SMBH in the universe (|Sesanall201ch . Tidal Disruptions of 
st ars by binary black holes have r e cently only been studied 
bv lChen. Liu, fc Magorriar] (|200Sl ); IChen et al.l (|2009| . |2010D 
and iLiu. Li. fc Chen (|2009l ) in the context of X-ray flares. 
They find, while the total amount of X-ray flares related to 
binary SMBH may be small (order 10%), they could show a 
special behavior in the form of bursts and interruptions of 
tidal disruptions. 

• The Bahcall- Wolf mechanism. 

In a collisional nucleus exchange of energy between stars 
drives the system towards an approximately steady-state 
distribution of stars around the SMBH in a relaxation time. 
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For a single stellar ma ss, this is f(E) ~ |_E| 1,/4 , p ~ r 7 1/4 
ijBahcall fc WohHl977T ). Since galaxies with collisional nu- 
clei probably always have M < 10 s M© the tidal disruption 
sphere is more relevant than the Schwarzschildradius. An- 
other condition for the Bahcall-Wolf solution is that rt is 
much smaller than r^m {\E t \ ^> GM./r^ m ), which is the 
case in real nuclei. It represents a 'zero-flux' solution. Nev- 
ertheless, in the numerical solutions, the steady-state flux is 
found to be small but non-zero. The flux is determined by 
the rate at which stars can diffuse into the disruption sphere 
at rt. 

The Bahcall-Wolf solution ha s been verified in a numbe r 
of other studies based on fluid (Amaro-Seoane et al.l 2004 ) 
or M onte-Carlo (|Marchant fc Shapirolll980l ; iFreitag fc Bend 
l2002h approximations to the Fokker-Planck equation. And 
it has been tested via direct NBody integrations, avoid- 
ing the approximations of the Fokker-Planck formali sm 
ijPreto. Merritt fc Spurzemll2004l ; iBaumgardt et al.|[20oi ). 

Observational measurements of the galactic center re- 
veal a stellar cusp, which appears to be flatter than ex- 
ected in a collisionally relaxed population around a SMBH 
Schodel et al"1l2007l ; |Po et al.|[2009l ; lMerritt|[20lb1 ') . Possible 
explanations for the absence of a cusp in the observed stars 
around Sgr A* include mass segregation (which leads to the 
formation of flatter cusps by lighter stars and steeper cusps 
by massive central stars, by relaxation), and the destruction 
of the e nvelopes of gian t stars in the densest parts of the 
cluster jPale et all 120091 ). It is not clear that Bahcall-Wolf 
cusps are present in any other galaxy however, both because 
relaxation times are generally 3> 10 10 yr, and also because 
they are difficult to be observationally resolved. 

Moreover, if binary SMBHs formed during galaxy merg- 
ers, they can destroy dense nuclei, as has been observed in 
the central d ensity profiles or 'mas s deficits' of bright ellipti- 
cal galaxies (|Merritt fc Szelll 12006). An important question 
is whether the existence of dense cusps at the centers of 
galaxies such as the Milky Way and M32 implies that no bi- 
nary SMBH was ever present, or whether a collisional cusp 
could have spontaneously regenerated after being destroyed. 

The relative importance of these and other various 
mechanisms, like physical collisions between stars, gas in- 
flow to the SMBH, and the nature of the dark matter that 
permeate galaxies and their possible interactions, are still 
not well understood. One of the most advanced physical de- 
scriptions including a detailed multi-phase treatment of in- 
terstellar matter and mechanical and thermal feedback due 
to stars and central b l ack h oles ( AGN) has been pr esented 
by ICiotti et al.l (|2009l , l2010f ) and IShin et al.l (|2010r ). Their 
approach, however, still lacks spatial and physical resolution 
compared to isolated galactic models. 



2 THE METHOD 

We define stellar accretion via the loss cone, which in a 
spherical galaxy is given by orbits, which specific energy 
(e) and angular momentum (J) lie within 

J 2 < Jl{e) = 2r t 2 [0(r t ) - e] ~ 2GM.r t (I) 

G is the gravitational constant, (f> the BH potential and 
Ji c denotes the loss cone size. Stars of mass m* and radius 



r* that come within a distance 

/2AO 1/3 



rt = a r* 



C2) 



of the central black hole with mass M, will be tidally dis- 
rupted and accreted. 100 % accretion efficiency is assumed. 
We include in this definition a free parameter a, for scal- 
ing of the tidal radius in our simulations (see description in 
Chapter 3). For M. < 10 8 M©, the tidal radius does not fall 
below the Schwarzsc hildradius, if we assume solar type stars 
(|Frank fc ReesHl976T ) 

Time scales considered here are of the order of the re- 
laxation time, or 



0.065 a 3 /(G 2 ?n*p In A) 



(■■>>) 



and In A is th e Coulomb logarithm ( Spitzcr 1987). We adopt 
A = 0.11A. (|Giersz fc Heggiej 1 1994^ 7 The relation between 
the dynamical or crossing time tdyn and t T is given by 



^dyn OC 



In A 
A" 



(4) 



where td yn — rja (r is a characteristic radius of the system, 
usually rhm). Here we point out the strong dependence of t T 
on N . We will specially scale our models in Sec. 3. 2. 



2.1 Influence and wandering radius 

Motion of stars surrounding the black hole in a certain re- 
gion are directly influenced by its gravitational field. This 
region is given by the influence radius of the BH, defined 
as the radius, where the mass in stars is of the order of M. 
(for an isothermal sp here M(< r n ) = 2M., Merritt 2003). 
iFrank fc Reel] (|l976T I define it by 



GM. 



(5) 



a is here the one-dimensional velocity dispersion. Though 
both definitions not always agree well, if the mass distribu- 
tion is known, the first one is easy to determine. 

The motion of a heavy particle in a sea of 
lighter particles can be treat e d as Brownian motion. As 
IChatteriee. Hernquist fc Loebl (|2002T ) point out, for King 
models with Wo ^ 3 (as used in the present study) equipar- 
tition of the black hole with its surroundin g core is a valid 
assumption |Merritt. Berczik fc Laurill2007l ). The wandering 
radius can be define as 



rwalk 



0.5 



1/2 



r,ib 



(6) 



where r n b is the NBody unit of length. This empirical rela- 
tion provides a decent fit for the measured black hole wan- 
dering during the whole simulated time in all considered 
models. This random walk might hinder the formation of a 
(7/4) cusp, as it stirs up the central region and smears out 
the eventual over-densities. Superposed on the random walk, 
it might have a non-vanishing mean velocity varying on 
much larger timescales. The black hole exerts high frequency 
oscillations, and has an additional dri ft, which is observed 
in th e center of mass vector as well ([Making fc Sugimotol 
Il987l ) . We correct this effect by subtracting every time step 
the center of mass from the black hole position. 
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2.2 Loss cone flux 

In a time t T , gravitational encounters between stars can ex- 
change orbital energy and angular momentum. Core collapse 
(shrinking of the core to zero size and infinite density) does 
not happen in galactic nuclei embedding a SMBH, since the 
central potential triggers core expansion and the central den- 
sity drops. 

The con c ept of a loss cone as introduced by 
iFrank fc Reel l| 19761 ) can be used to identify stars on an 
orbit penetrating the black holes Roche lobe. As these loss 
cone orbits are disrupted within an orbital period, the stel- 
lar mass supply must run out immediately if they are not 
replenished by relaxation. The loss cone can be defined as an 
angle like variable giving the half aperture of the black hole 
as seen from the stars distance. Trajectories with a smaller 
aperture have a peribothron r t and will be lost. 

In a potential dominated by the black hole, taking ad- 
vantage of the Keplerian velocity profile, 



GM. 



(7) 



and using Eq. [T] for e < e t l|Frank fc Reeslll976h . the loss 
cone angle 6\ c — v\ c /<J evaluates to 



1/2 



(8) 



In a steady state spherical system, the only process re- 
filling the loss cone is scattering of stars due to distant grav- 
itational encounters. For a general r~ v density cusp and a 
given by Eq. [7] the relaxation time within the black hole 
sphere of influence becomes 



tr = 0.338- 



M. 



3/2 



„>7-3/2 



'GWmtPnlnArZ' ^ 

where r n and p n are radius and mass density at the influence 
radius, and the angular diffusion per orbital period 9o = 
(tdyix/tr) 1 turns out 



I = 2.960 r 3 ~ v r n , M~ z m i ,p n In A. 



(10) 



Comparing the two angles, one customarily defines two 
regimes: The empty loss cone or diffusive regime for small 
r where 9\ c > 8u and the full loss cone or pinhole regime, 
where stars can move t hrough the loss cone within one or- 
bital period (#i c < 8n). iLightman fc Shapiro] (| 19771 ) showed 
that the flux into the loss cone peaks at 9u = 0\ c , which 
defines the critical radius 



0.225 Mt 



ft 



m* p n In A 



(11) 



for the assumption of a power law density profile within 
the black holes radius of influence. Hence the last relation is 
valid for r C rit i$ Him, whic h holds for most of the 5 1 elliptical 
gala xies in the sample of jWang fc Merrittl (|2004l ). 

ILightman fc Shapiro! (1 19771 ) showed that the integrated 
number flux per orbital period, can be approximated to 
match its value at the critical energy, 



F(E) ~ F E (E crit )\E crit 



(12) 



where Fe gives the flux of stars at E — E CI i t . A simplified 
expression of the expected disruption rate for the steady 
state solution (77=1. 75), can be obtained by applying the 



scaling of N(r) tx r 5//4 , and t r oc a A /n(r) oc r 17//8 in Fe oc 
N(r)/t(r), and E oc r" 1 , an d the scaling of r(t) oc t 2/3 
during self similar expansion (|Shapird Il977l ). It results in 

We use an equivalent expression jFrank fc Rees|[l976r i 



Acrit OC 



4nr 3 6f c (r)n(r) 



3 tdyn(r) 

for a n = 7/4 power law cusp and tdyv 



r 3/2 



(13) 



to derive 



A 7/4 oc 6.39VGlnA B / 9 M- n ' la rt /9 ml 0/9 nl 4/ \t 9 /W (14) 
And in physical units, 



A7/4 oc 



6.89 x 10" 



Myr 

( M. 

V 1000M e 



In 



-25/24 



M. 

2ml 



no 
pc _: 



5/9 



Rq 



4/9 



Mo 



14/9 



49/18 



26/27 



(15) 



This equations will help us to properly scale our simu- 
lations, as described in the next section. 



3 SIMULATIONS AND RESULTS 

Evolution of dense stellar systems harboring growing black 
holes is studied using direct NBody methods with an imple- 
mented tidal disruption procedure. The black hole is treated 
as a heavy particle with an initial mass of M./Mtot = 0.01. 
Particle numbers are N = 10 000 to N= 100 000. Simula- 
tions were run in parallel with up to 128 processors on the 
RZG Power 6 Machine in Garching (Munich), which is a 
facility of the DEISA project; the 40 nodes Kolob CPU- 
Cluster (Mannheim Germany), the 85 nodes Laohu cluster 
(Beijing, China) and the GRAPE cluster Titan at ARI-ZAH 
(Heidelberg, Germany), which has been recently upgraded 
with GPU cores. In our implementa tion, by using the MPI 
parallel NBo dv6++ JSpurzeml[l999l ) and the parallel GPU 
code y;GPU (|Harfst et al.ll2007 ). we let the star proceed to 
its peribothron where the accretion then takes place. In the 
NBody6+- 1- code, we tak e advantage of the neighbor scheme 
|Ahmad fc Cohenl Il973l ). and find candidates for the dis- 
ruption by only searching the black holes neighbors - thus 
reducing the computation al overhead. These simulations in- 
corporate zero softening |Harfst et al.l 120071 ) and have an 
energy conservation AE / Eq < 10 -4 , even after 10 4 NBody 
time units. For the high N (up to 100K) simulations we use 
the parallel ipGPU code, ready for use with GPU clusters 
which includes a softening parameter of 10 -5 . The energy 
conservation is of the same order as in the NBody6H — h code. 

In order to challenge the analytic theory of loss cone 
diffusion by direct NBody simulations, we have to ensure 
that the simulations probe the dynamics of interest for real 
galaxies - the empty loss cone regime. It is convenient to 
identify the physics at work by its dominating time scale, 
which in our case demands a separation of the loss cone 
depletion time tout = £d yn from the loss cone refilling time 
t- m — 0f c t T . Hence the condition to have an empty loss cone 
is 



rt 

^dvn tl *C 

r 



(16) 



Evolution of growing black holes in axisymmetric galaxy cores 5 



For a in-depth d iscussion of tj n we like to r efer to the gas 
model studies bv lAmaro-Seoane et alj (|2004j. A straightfor- 
ward way to satisfy the above relation is by increasing the 
particle number N following Eq. [4] However, as the 0(N 2 ) 
scaling of direct algorithms transforms to 0(N 3 ) for relax- 
ation processes, it was a challenge already for our N=100K 
runs, specially because of the long integration time required 
for our purposes. Additionally, we highlight that the inclu- 
sion of a heavy black hole particle leads to a wide-spread 
time step distribution with few very short stepped particles 
close to the black hole. This hampers scalability and signif- 
icantly increases the over-all integration time compared to 
the standard benchmark cases. Thus, it is still very difficult 
to obtain models of N ~ 10 6 on todays general purpose high 
performance computers, but we expect in a soon future to 
be able to perform such runs in our new GPU clusters. 

Given the limitations on the particle number and to fur- 
ther separate the time scales, we introduce the magnification 
factor a (Eq. [2| to vary the tidal radius in our NBody sim- 
ulations with this free parameter and in order to determine 
the scaling behavior of the results as a function of a, which 
ranges between 1 to 1000 (or ~ 10~ 5 rh m < r t <~ 10~ 3 rh m ). 

It allows us to improve the performance and as well a 
deeper study of time dependent stellar accretion. Our model 
NBody units (jHeggie fc Mathieulll986h are scaled to G = 1, 
r^b = 1 P c , and m, = IMq. With this scaling the stellar 
radius for the single mass runs turns out to be r, = 2.52 x 
10 _8 r n b the tidal radius then turns out as r t (0) = 2.52 x 
10~ 7 r n b for a mass ratio of 1000:1. We are neglecting in this 
study of single-mass systems, the stellar evolution during 
the whole simulated time. We are presenting evolution of 
multi-mass systems in a forthcoming publication, where we 
will fully incorporate the stellar evolution time scale, which 
becomes short enough to play a significant role. 

The black hole is located initially at the center with zero 
velocities and the mass of tidally disrupted stars is added 
completely to its mass every accretion event. In order to 
provide an actual sink in phase space to drive the diffusion, 
particles need to be removed from the simulation when en- 
tering the tidal radius. The subsequent accretion is bluntly 
modeled as a perfect inelastic collision with the black hole 
particle, where the star is fully accreted and linear momen- 
tum is conserved. The "equations of motion" read 

M' m = M, + m* (17) 

r'. = J^(M.r.+m^) (18) 

v'. = M^-J^+m^). (19) 

As initial galaxy mo dels, we use rotating King Models 
|Einsel fc SpurzerrJ 119991 . where we added the black hole 
particle of mass M./M to t = 0.01 in the center. For com- 
pleteness we employ Wo = (3, 6) and uo = (0.0, 0.6, 0.9) ax- 
isymmetric King models. Wo — 3 King models have larger 
and denser cores than Wo = 6 King models. In the text, we 
refer to non-axisymmetric models (ojo = 0.0) as an approach 
to spherically symmetric systems, since the first don't have 
any flattening due to rotation. We use, complementary to 
our NBody models, Fokker-Planck approximations to our 
problem, with a seed BH M. = 10~ 5 and N = 10 s to 
properly define Ji c , as in Eq. Q] Here is the stellar radius 
r* = 2.52 x 10~ 8 r c , and thus r t (0) = 2.52 x 10~ 7 r, where r 



Table 1 . Overview of the performed NBody runs 
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Col. (1) is an identifier for the run, Col. (2) is the King parameter 
of the initial model Wo, Col (3) is the rotation- parameter in the 
King-model i^iq, Col. (4) is the initial number of particles, and 
Col. (5) is tidal radius magnification factor a 
Notes: computer cluster where simulations were performed, 
(a) NBody6++ in RZG, (b) ^GPU in Titan, (c) i^GPU in Kolob, 
(d) i^GPU in Laohu (see text for a brief description of the hard- 
ware) 



is the scaling radius in our King Models, equal to the core 
radius. The value of the initial rt corresponds to a mass ra- 
tio M./m* — 1000. Table [T] summarizes the parameters of 
the performed NBody and Fokker-Planck runs. 
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Figure 1. Evolution of the Lagrangian radii for axisymmetric King models 16KR3c (loft) and 16KR6c(right). Time is given in units of 
i r h(0) and radii in units of r^ m (0). The system expands after t cc due to the presence of the BH. Evolution of r n (green line), r cr ; t (yellow 
line), and r^ik (light blue line) is shown. The initial M./Aftot = 0.01 



3.1 Cluster evolution 

The negative heat capacity of self gravitating systems leads 
to the pr ocess of core c ollap se. In the stand ard picture pro- 
posed by IHenonl (|l965l ) and lAarsethl (jl973h . the singularity 
is avoided due to the formation and subsequent hardening 
of tight central binaries. The presence of a star-disrupting 
black hole can act as an energy source just like binary hard- 
ening. On the event of star-disruption, a bound object is 
removed from the system, which looses (negative) binding- 
energy, (energy conservation is established, if the inner en- 
ergy of the black hole is considered) the refore the system 
gains energy and expands. IShapirol (jl977t ) studied the com- 
bined effect of star disruption and the change of density in a 
homological model. He found a self-similar expansion of the 
core-radius (r c ) according to 



r c (t) 



cx [1 + g{M.) t] 



2/3 



(20) 



r c (t cc) 

where t cc is the collapse timeQ, and g{M m ) also depends on 
the initial- and on the "minimum-" core radius and respec- 
tive density, as well as on the relaxation time. Given the 
similar underlying physics it is unsurprising that the time- 
dependence is of the same ty pe (oc t 2 / 3 ) a s in the binary 
harde ning case considered by IHenonl (|l965l ) and iGoodmanl 
d 19841). 



We show the Lagrange radii (radii containing the given 
percentage of the initial total stellar mass) for the runs 
16KR3c and 16KR6c in Figure [T] We are using a small par- 
ticle number, what permits us to obtain longer evolutionary 
times, in relaxation time scales. 

After a less pronou nced collapse co mpared to the 
case without black hole dKim et all 120081 ). the La grange 
radii expand self similarly according to the r oc t 2//3 
law. This behavior, which has been also seen in gas mod- 



the collapse time corresponds to the time at which central den- 
sity grows to infinity, and core radius shrinks to zero. In systems 
with an energy source, like here, it is the time at which the con- 
traction phase is halted and reversed 



els (|Amaro-Seoane et al.l |2004|) and axi symmetric Fokker- 
Planck models ( Fiestas fc SpurzemlbOlOl ) , is here verified in 
a self consistent direct NBody simulation for axisymmetric 
systems. 

As we can see, Lagrangian radii give us a qualitative 
description of the interaction of a growing BH and the stel- 
lar mass shells. Initially, Lagrangian radii are dominated by 
core contraction and the BH mass growth is slow due to the 
low central density. Later, density grows due to gravitational 
instabilities (the Lagrangian radii shrink) and the collapse is 
halted and reversed (mass shells re-expand), while the grow- 
ing BH potential, dominates the system. The axisymmetric 
models of Fig[T]show a similar evolution as the well studied 
spherically symmetric case. 

In Fig[T]we observe that t cc is shorter than the expected 
value for systems without BH, which are £ C c/t r h ~ 16 for 
Wo = 3 models, and t cc /£rh ~ 10 for Wo = 6 models. This 
effect is induced by the magnification of the tidal radius and 
the high initial M., which enhance stellar accretion, acceler- 
ating the reverse of collapse and further expansion. Nonethe- 
less, it does not affect the obtained self-similar expansion. 
Note that the influence radius approaches the initial rh m 
during the post-collapse phase. r cr it (Eq. Ill[) is also shown 
in this figures, as well as r wa ik (Eq. [6]), which decreases with 
increasing M, , and becomes more than an order of magni- 
tude smaller than r cr it , and more than 2 orders of magnitude 
smaller than r n . 

Spherically symmetric systems with BH are known to 
develop steady-state solutions in relaxation time scales. 
These solutions have characteristic density and velocity dis- 
persion profiles. The density distribution scales theoretically 
with radius as r~ 7//4 and the velocity dispersion as r^ 1 ^ 2 . 
We follow the evolution of our models in this time scale and 
compare them with the standard theory. 

Fig. [5] shows the evolution of the density profile for the 
models lOOKRlc, 100KR4c (non- rotating), and 100KR3c, 
100KR6c (rotating). We use here the highest particle num- 
ber in order to obtain a better resolution of the cusp. Note 
that the cores are initially flat and the time needed to build 
the cusp is of the order of t T h- The influence radius (green di- 
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Figure 2. Evolution of the density profile of non-rotating models lOOKRfc and 100KR4c (top), and rotating models 100KR3c and 
100KR6 (bottom). Symbols arc: Influence radius o and critical radius A. In post-collapse, the radii are well separated and the profile 
approach the expected zero flux solution between r rcr ;t and v n . 



amonds) moves outwards while BH mass grows, and the crit- 
ical radius (orange triangles) appears in our resolution range 
before the cusp is formed. r wa ik is smaller than the minimum 
radius in this figures. In post-collapse, these radii are well 
separated and the profile approaches the expected zero flux 
solution between r cr i t and r n . Note that systems with larger 
cores (Wo = 3.0) evolve slower than Wo ~ 6.0 models, reach- 
ing the latter their final density cusps in shorter times. 



During post-collapse the system expands and the cut- 
off radius in our models extends up to ~ 10 rh m - Closer 
to the center, one can see that at later times, the critical 
radius grows and the cusp is shallower inside this radius, 
due to the effective stellar accretion around the BH and the 
growing tidal radius, which perturbs the formation of a cusp 
at the very center. Although the critical radius is resolved 
in our simulations, it contains only few dozens of stars. We 
are performing higher N (Js 256K) simulations, which will 
provide a more accurate measurement of the central cusp 
and specially of the empty loss cone region. 



3.2 Disruption rates 

In our models, stellar accretion is driven by small angle, two- 
body encounters, which under the influence of the BH grav- 
itational potential, causes that some stars lose energy and 
move closer to the BH being eventually consumed. During 
the contraction phase, high central densities are expected to 
trigger higher BH mass growth rates. Thus, maximal rates 
occur close to t cc , when angular momentum and energy dif- 
fusion are most effective. 

As previously discussed, in our NBody models, stars in 
orbits of J < Jic, which reach their apocenters inside the 
tidal radius define an accretion event. In our axisymmet- 
ric Fokke r-Planck models, stars in orbits of J z < J Zj i c are 
accreted (jFiestas fc Spurzemll201ol ). being only energy and 
J z conserved quantities. Since the other angular momentum 
components are not conserved, accretion can be artificially 
enhanced in this models, specially during the initial evo- 
lution, before self-similar expansion sets on. Regarding ini- 
tial conditions, our Fokker-Planck models use comparatively 
smaller BH seeds (^21 ^ icr 5 ), while we fix this value in 
the NBody models to 0.01. 

Since we expect self-similar evolution during core ex- 
pansion, independent on initial conditions, as shown in 
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Figure 3. Left: Non-scaled disruption rates for NBody (16KRla,b,c) and Fokker-Planck models (FPKR1). The influence of the magnified 
rt in the evolution is clearly seen. The solutions approach the Fokker-Planck results for a = 1. Right: maximal disruption rates vs. 
parameter a for all models. Filled symbols are Wo = 6 models, empty symbols correspond to Wo = 3 models. The results converge to 
the ideal case (o = 1) as in the Fokker-Planck approximation. 
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Figure 4. Disruption rates in units of fractional mass per relaxation time for King models 16KRlc, 16KR4 (non-rotating, top) and 
16KR3, 16KR6c (rotating, bottom). Nbody results are black lines, green lines are Fokker-Planck models and the expected time dependence 
for a Bahcall-Wolf cusp (oc t~ 1 - 25 ) is shown by the red line. The predicted rate following Eg. 1 141 is shown as a blue line. The rate peaks 
after a few t T h and decreases afterwards. The models follow the analytical predictions in the post-collapse phase. 
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Figs. [T] and [2] at this stage both methods should be 
comparable, if radial and time units are properly re- 
scaled. We check this by scaling rates, correcting the 
changes induced in Eq. [14] by M,, m* = 1/N and r t in 
each model. First we bring rates to common dimension- 
less units ((AN/N)/(At/t T h)), and we calculate the factor 
JVmodeiA/^VmodciB, by replacing the correspondent values of 
Mm, rrii, and rt of each model. r t enhances the disruption 
rates by a factor of ~ 1 to 22, (for a = 1 to 1000). A vari- 
ation of N between 10K and 100K modifies the stellar mass 
rrii, as l-6x 10" 5 , and M./m* changes like ~ 100-1000. A 
(-7/4) cusp in the density profile is assumed. We scale finally 
the results obtained by model B to model A by multiplying 
the first with the previous factor. 

On the other side, as previously discussed, the larger 
seed BH in the NBody models triggers core heating and 
enhances accretion events, with the consequence of an ear- 
lier bounce of core density. It modifies the time at which 
density reaches its maximum, and maximal disruption rates 
appear. We can observe this effect in the evolution of La- 
grangian radii (Fig. [TJ. Thus, self-similar expansion during 
post-collapse can be compared by bringing the collapse times 
of the compared models together. 

Fig. [3]left shows disruption rates before scaling. Note 
that peak rates are higher for higher a and NBody models 
approach Fokker-Planck models by decreasing this param- 
eter (the 'real' tidal radius is given by a — 1). Since max- 
imal disruption rates dominate the stellar accretion events 
over time, they can give us the extent of the influence of 
the different initial parameters and kinematics (rotation) 
used in our models. Figure [3]right shows maximal disrup- 
tion rates (7V max ) obtained for all models in dimensionless 
units (dM / Mtot) / (dt /t T h) against the parameter a, as used 
in the simulations. This plot gives only approximately the 
7V(a)-dependence, since the m*-, and M.-dependence are 
also present. Additionally, the different peak rates at con- 
stant a, are a consequence of the better resolution reached 
in high-N Nbody runs, which show higher peaks. Moreover, 
iV max , for smaller values of a converges to rates correspond- 
ing to a = 1 (actual r t ), as we show by plotting the maximal 
rates of our Fokker-Planck models (black dots in Figure [3] 
right). The measured values of jV max are shown in Table [2] 
Mainly models using a — 1000 and ljo = 0.0, 0.9 were here 
selected for comparison. Table [2] shows as well iV^J«j scaled 
to a = 1, as in our Fokker-Planck models. We correct for 
the N oc r\ ^-dependence, as in Eq. 1141 

Fig. [4] shows the disruption rates for models of King 
Parameter Wo = 3.0 (left) and Wo = 6.0 (right). In 
each figure Nbody results (black lines) together with scaled 
Fokker-Planck approximations (green lines) are plotted. Ad- 
ditionally the expected time dependence of iV oc t -1,25 , for 
r\ — 1.75, following Eg. 1121 is shown (red line). We introduce 
as well the expected rate in our simulations for a (7/4)- 
cusp by using Eq. [14] (blue line). Rates are averaged over 15 
events and smoothed over 20 t n b. Rates initially increase fast 
and reach their maxima close to t cc , decreasing afterwards 
during the expansion of the system. From these results, we 
can conclude that in the post-collapse phase, our NBody 
models follow the predicted evolution and agree with the 
Fokker-Planck results during self-similar expansion. We can 
further test both models by comparing scaled maximal dis- 
ruption rates, as shown in column 5 of Table [2] In physical 



units, for a galaxy core of M ~ 10 9 A/q and a relaxation 
time of 10 Gyr, we obtain N max = 1.3 ±0.2 x 10 -4 M o yT -1 . 
Here we average the results of our (Wo = 0.6, ojo = 0.0)- 
models and correct the r t (a and N)-dependence as in Eq. 1141 
and Eq. [2] The corresponding Fokker-Planck models give 
Amax = 1.1 x lO _4 M0yr _1 , in agreement with the NBody re- 
sults. Thus, we can conclude that we obtaine a non-trivial re- 
sult through all the post-collapse, where NBody and Fokker- 
Planck models agree with each other, which is of particulary 
importance, since both methods use an equivalent but not 
identical treatment of accretion and because definition of 
relaxation time (our Eq. [St in galactic nu clei can be eventu- 
ally over-estimated (|Madigan. et al.ll201ll ). This means that 
we can from these data make predictions for the star accre- 
tion rates and other parameters of dry galactic nuclei, which 
are independent of the previous history. And it shows, that 
relaxation processes might play an important role in the 
growth of SMBH. 

By comparing non-rotating with rotating models in Ta- 
ble [2] column 4, we can observe that the latter show slightly 
higher peak rates, for constant N. Maximal rates vs the rota- 
tional parameter u>o for all models are shown in Fig. [5] top. 
Low N runs show roughly constant peaks, while the better 
spatial resolution reached by high N runs, specially at r C rit 
(as seen in Fig. [2]| , show higher disruption rates for higher 
cjo- The higher peaks reached by rotating systems lead to 
the BH masses shown in Table [2] Final BH masses (M m j), 
are measured at times, which are kept the same for constant 
N and Wo to facilitate the comparison between the models. 
We observe specially in the N=100K runs, that non-rotating 
models reach a similar final mass, since rotating models con- 
verge to a ~ 20% higher mass, independent of Wo - It implies 
a direct influence of rotation in the process of stellar accre- 
tion. 

In order to understand how this effect is triggered, 
Fig. [5]middle shows the distribution of J x , J y , J z of accreted 
stars for the model 100KR6 (rotating). Values are fractions 
of the total number of accreted stars. All distributions peak 
at J = 0, as expected. The contribution of J x , J y to ac- 
cretion is symmetric with respect to J = 0, while a higher 
fraction of prograde rotating stars (J z > 0) is consumed by 
the BH. This was expected, since our initial models rotate 
in the positive direction of J z . Fig. [5]bottom shows the dis- 
tribution of J z of accreted stars for the model 100KR6 (ro- 
tating) in comparison to the non-rotating model (100KR4). 
We observe that rotating models show an excess of accreted 
prograde rotating stars. 

We show in Fig. [6] the radial distribution of the semi- 
major axis, normalized by r cr i t , of accreted stars in models 
100KR4c (non-rotating) and 100KR6c (rotating). A larger 
fraction of accreted stars in the rotating models have semi- 
major axis few times larger than rcrit- Distribution of ac- 
creted stars in the non-rotating mode l peaks at a/r c rit ~ 1 
as ex pected from the loss-cone theory (jLightman fc Shapiro! 
Il977t ). The distribution of eccentricities is similar in both 
models (Fig. [6] bottom). It shows a steep maximum above 
e = 0.995. From Figs. [5] and [6] we can conclude, that the ex- 
cess of accreted rotating stars, origins mainly from regions 
outside r cr i t . According to loss cone theory, diffusion of en- 
ergy and angular momentum is higher than the loss cone 
size in this region (full loss cone regime), and orbits in this 
region are able to escape the loss cone in a dynamical time. 
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Figure 5. Top: Maximal disruption rates vs. rotational parame- 
ter too for all models. Symbols are like in Fig. [3] right. Specially 
high N runs show how rotation influences the disruption rates. 
We keep a constant parameter a = 1000 for comparison. Mid- 
dle: Distribution of J x , J y , J z of accreted stars for the rotating 
model 100KR6. A higher fraction of progradc rotating stars is ac- 
creted. Bottom: Distribution of J z of accreted stars for the mod- 
els 100KR4 (non-rotating) and 100KR6 (rotating). The excess 
of accreted prograde rotating stars leads to the obtained higher 
masses. The maximum J z j c , for the largest energy of the stel- 
lar orbits, is indicated in the middle and bottom panels with a 
vertical dashed line. 



Figure 6. Top: Histogram of semimajor axis over critical radius 
of orbits of accreted stars in models 100KR4c (non-rotating) and 
100KR6c (rotating). Most accreted stars in the rotating models 
have semi-major axis few times larger than r cr it- Distribution of 
accreted stars in the non-rotating model peaks at a/r cr ; t = 1 as 
expected. Only bound orbits are shown in this figure. Bottom: 
Histogram of eccentricities for the same accreted stars. Distribu- 
tion in both models is similar. A contribution of unbound orbits 
is present. 



In order to obtain a larger r CI n, as in our rotating models, 
tin ~ tfLtr should be larger, or £ out ~ 8^t r shorter in these 
region Q This hints to a breakdown of the classical theory 
which depends on the conservation of angular momentum 
that is not given in the rotating models. An enhanced re- 
laxation in a system that is supported by bulk rotation, as 
compared to a non-rotati ng one which is pressure supported, 
could explain this effect. (jAthanassoula. Vozikis fc Lambert! 
l200ll) . A decrease in the velocity dispersion for r < r n , which 
could modify substantially t T oc a 3 within the BH influence 
zone, was nevertheless not detected in our models. It is thus 
tempting to interpret this behavior in terms of the additional 
presence of orbits with non-conserved J x , J y angular momen- 
tum (e.g. box orbits). To further investigate this effect, an 



2 the crossing point in Fig. 13 of lAmaro-Seoane et al. is 
shifted to the right in our rotating models, wenn compared to 
non-rotating models 
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Table 2. Results of the simulations 



Identity 


W 


[JO 


jymax 


7V ma , x . 

scaled 


M ;f 


*/*rh,f 


16KRlc 


3.0 


0.0 


.0252 


.0008 


0.2804 


31.5 


16KR3c 


3.0 


0.9 


.0254 


.0008 


0.2788 


31.5 


16KR4c 


6.0 


0.0 


.0193 


.0008 


0.2836 


31.5 


16KR5c 


6.0 


0.6 


.0231 


.0010 


0.2803 


31.5 


16KR6c 


6.0 


0.9 


.0232 


.0010 


0.2785 


31.5 


32KRlc 


3.0 


0.0 


.0323 


.0012 


0.2048 


10.2 


32KR3c 


3.0 


0.9 


.0358 


.0012 


0.2102 


10.2 


32KR4c 


6.0 


0.0 


.0358 


.0013 


0.2624 


13.0 


32KR5c 


6.0 


0.6 


.0359 


.0013 


0.2776 


13.0 


32KR6c 


6.0 


0.9 


.0358 


.0013 


0.2799 


13.0 


64KRlc 


3.0 


0.0 


.0535 


.0018 


0.1592 


5.5 


64KR3c 


3.0 


0.9 


.0468 


.0016 


0.1791 


5.5 


64KR4c 


6.0 


0.0 


.0528 


.0017 


0.2283 


6.6 


64KR6c 


6.0 


0.9 


.0601 


.0020 


0.2647 


6.6 


lOOKRlc 


3.0 


0.0 


.0600 


.0018 


0.2478 


6.9 


100KR3c 


3.0 


0.9 


.0800 


.0020 


0.2785 


6.9 


100KR4c 


6.0 


0.0 


.0611 


.0019 


0.2542 


6.9 


100KR6c 


6.0 


0.9 


.0850 


.0026 


0.3175 


6.9 


FPKR1 


3.0 


0.0 


.0008 


.0008 


0.0011 


30 


FPKR3 


3.0 
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Col. (1) is the identifier for the run, Col. (2) the King parameter 
of the initial model Wo, Col (3) the rotation-parameter in the 
King-model ujq , Col (4) are the maximal disruption rates in units 
of {dMm /Mtat)/(dt/tih) i Col (5) are the rates scaled to a = 1, 
Col. (6) is the final black hole mass, and Col (7) is the time in 
relaxation time units at which the final mass was measured. 

orbit study of the accreted stars should be performed, which 
we must leave open for future research. 

3.3 Rotational velocity 

In Fig [7]top we show the initial distribution of the veloc- 
ity component correspondent to J z , of stars for the model 
100KR6c. The initial rotating King Models show a maxi- 
mum of rotation (u ro t,max) at around rhm and central rigid 
body rotation. We plot velocities of bins of 5 stars, in or- 
der to get a more detailed stellar distribution, and average 
over 5 NBody time units. Vertical dashed lines mark the 
radius of influence of the BH and r cr it in units of r^ m at 
the given time, as indicated. The rotation profile from bins 
of 50 stars is overplotted (red dots). The evolved distribu- 
tion of rotational velocities for the system after relaxation 
(~ 3 t r h) is shown in Fig. middle. As a consequence of an- 
gular momentum transport through gravitational scatter- 
ings, the original i> ro t,max decreases during the evolution. 
In the model w r ot,max(icc)/wrot,max(0) ~ 0.85. Additionally, 
one can observe that u r ot,max moves inwards, in a region 
close and inside the influence radius, leading to an increas- 
ing of rotation in this region, with respect to the initial 
configurations. With constant ui = 0.9, as an initial rota- 
tional parameter, this effect is more pronounced in concen- 
trated models (100KR6c, Wo = 0.6) than in models with 



larger cores (100KR3c, Wo = 0.3), because the first con- 
tain initially more amount of rotational energy with respect 
of kinetic energy (T ro t/2kin ~ 0.3) than the second ones 
(Trot/Tkin ~ 0.1). After the system relaxes, central stars ro- 
tate with velocities in average larger than the original max- 
imum, initially located in the outer parts of the system. 

Consider the region r cr it < r < r n , where stars are 
dominated by the BH central potential. Stars populate al- 
ways more this region in time (Fig. (7]middle), while they 
interact dynamically with other stars. A fraction of them 
will be disrupted, when their orbits reach r < r t . As we 
discussed, the excess of disrupted stars in rotating models 
comes mainly from this region, and is dominated by stars 
in orbits with positive J z . It is not surprising, since the con- 
centration of rotating stars in the BH zone of influence is 
triggered by the initial configuration in our axysimmetric 
systems. Moreover, the growing central density, caused by 
gravi tational and gravogyro instabilities (|Einsel fc Spurzeml 
Il999h together with angular momentum transport, enhances 
this effect, specially before expansion sets on. 

The few counter-rotating stars observed in Figure [7J 
middle, are remaining stars from the initial distribution, 
which together with pro-rotating stars lead to no rotation 
in the very center. The enhancement of rotation inside the 
influence radius builds a wide maximum of rotating stars in- 
side r n , which is now close to rh m - At a radius inside ~ r cr it 
one finds only a few dozens of stars (in our N=100K runs), 
and is difficult to talk about rotation in this region. More- 
over, in the region r cr i t < r < r n a rotating core can be 
detected. By comparison with Fig. [2] we observe that this 
is the region where the cusp forms, and the BH potential 
dominates the stellar environment. 

In order to investigate how strong rotation dominates 
the dynamics at this evolutionary stage, we show the rate 
of rotational velocity vs. one dimensional velocity dispersion 
(urot/c) for the same relaxed axisymmetric model (Fig. [7] 
bottom) correspondent to the previous figures. This param- 
eter shows the relative importance of rotational vs. pressure 
supported kinematics as used in observational studies of el- 
lipticals and galaxy bulges of spirals. Here we used velocity 
bins of 50 stars and averaged over 50 time steps in order to 
get a better defined profile. In relaxed systems, the initial 
peak at ~ r^m still dominates and becomes wider. Although 
the profile continuously decreases towards the center, an en- 
hanced Wrot/cr inside r n with respect to the initial configu- 
ration can be observed. It is specially interesting since the 
BH gravitational potential builds a velocity dispersion cusp 
(or. r -1 / 2 ), which requires a strong increase of v rot in order 
to be detected. 



4 CONCLUSIONS 

Lo ss cone theory as d eveloped in the cla s sic p apers 
of iFrank fc Rees] (1 1976)) ; iLightman fc Shapiro! (| 19771 ) and 
ICohn fc Kulsrudl (1 19781 ) can be used to estimate feeding 
rates for SMBHs in galactic nuclei. Nonetheless, total con- 
sumption occurs from orbits that could extend beyond the 
black hole influence radius, hence the contribution of the 
stars to the gravitational potential cannot be ignored. These 
orbits interact with central stars and are able to interchange 
energy and angular momentum in relaxation time scales. In 
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Figure 7. Top: Initial distribution of rotational velocity of stars 
in the rotating model 100KR6c. Rotational velocities are aver- 
aged in bins of 5 stars and over 5 NBody time units. Vertical 
dashed lines mark r cr it and r n in units of r^m at the given time. 
The averaged profile is overplotted (red dots) and error bars are 
included. Middle: rotational velocity of stars after relaxation. A 
rotating core can be detected for r cr jt < r < r n . Bottom: Initial 
and final rotational velocity over velocity dispersion for the same 
model. Growth of v TO t/cr inside r n is less pronounced because of 
the central cusp of a. 



this work, we investigate accretion rates in axisymmetric 
systems by using direct NBody and Fokker-Planck simula- 
tions, harboring a star accreting growing black hole, which 
evolves dynamically in time scales of relaxation. 
Our main results are: 

• The stellar distribution is strong influenced by the in- 
terplay between diffusion of energy and angular momentum 
(gravitational instabilities) and the feeding of the black hole 
influence zone by stars in high energy, low J (eccentric) or- 
bits. These systems undergo core collapse (reach a central 
density maximum at t cc ) in the presence of a star-accreting 
black hole. The growing central BH potential dominates al- 
ways larger zones of the system, with a growing influence 
radius, which reaches almost the half-mass radius during 
the post-collapse phase. 

Axisymmetric (rotating) systems, like the spherically 
symmetric, reach in our simulations steady-state solutions 
during the post-collapse phase. Systems with smaller cores 
(Wo = 6.0) reach self-similar expansion in shorter times than 
systems with larger cores (Wo = 3.0). For relaxed systems, 
the BH r wa ik is about one order of magnitude smaller than 
r cr i t , which is well resolved in our simulations, and itself 
about two orders of magnitude smaller than r n . Thus, these 
systems fulfill the conditions for the loss cone theory. 

Furthermore, we are currently investigating the conse- 
quences for relaxation times and specially for the loss cone 
theory in multi-mass models. Mass segregation time scales in 
real systems are short enough to lead to a faster relaxation in 
the central parts, where the more massive stars concentrate, 
and consequently, to an earlier formation of cusps by these 
central stars (a work to be presented in a forthcoming publi- 
cation). This is specially important, since galactic nuclei are 
often less than one relaxation time old, and stellar density 
near the SMBH seems not to have the Bahcall-Wolf form, 
although observations need to resolve the influence radius to 
be able to detect them. Additionally, relaxation times can 
be themselves much shorter, the smaller the radius of the 
system, and the higher the central densities. 

• We measured tidal disruption rates in axisymmetric 
NBody models and compare them to 2D Fokker-Planck re- 
alizations. Tidal star accretion acts as an indirect heating 
source reversing core collapse, like in isolated star clusters it 
is the role of hard binaries. In our low N, NBody runs with 
high initial seed black hole masses the indirect heating is 
stronger, so the density maximum, and thus the peak rates 
are wide. For high N or small seed black hole masses the 
heating is small, a very peaked central density and disrup- 
tion rates occur (like in the FP models), as shown in Fig. [4] 

In order to investigate self-similar evolution during the 
post-collapse phase, we apply a scaling procedure to our 
models, which use partly different initial conditions and 
we obtain a non-trivial result for the rates through all the 
post-collapse phase where NBody and Fokker-Planck mod- 
els agree with each other. This means that in this phase we 
can from these data make predictions for stellar disruption 
rates and other kinematical parameters of dry galactic nu- 
clei, whic h are independent o f the p revious history (as also 
shown in iFiestas fc Spurzeml (|2010h 'l. even independent of 
whether the black hole has grown earlier by gas or star ac- 
cretion. A statement, which is confirmed here by using direct 
NBody realizations. 



Evolution of growing black holes in axisymmetric galaxy cores 13 



We found that disruption rates and BH masses are in- 
fluenced by axisymmetry/rotation, in the way that rotation 
leads to higher peak rates and higher M.j. The excess of 
accreted stars, origins mainly from prograde rotating stars, 
located in regions outside r rit- This hints to a breakdown 
of the classical theory, given by the rotating models, which 
could be interpreted as a consequence of the presence of box 
orbits with non-conserved J x ,J y angular momentum. To fur- 
ther investigate this effect, an orbit study of the accreted 
stars is necessary. This is of importance, since galactic nu- 
clei need not be even axisymmetric. In a triaxial nucleus 
containing centrophilic orbits, the mass in stars on orbits 
that intersect the SMBH's capture sphere can be enormous, 
much greater than M,, so that the loss cone is never fully 
depleted. Galactic nuclei sometimes undergo catastrophic 
changes, due to galaxy mergers, in-fall of star clusters or 
black holes, star formation, etc, all of which can substan- 
tially affect the feeding rate on both the short and long 
terms. 

We apply, for illustration, disruption rates given by Eg. 1151 
to the galactic center, by using M. = 3.3 x 10 6 , r n = 1.65, 
po = 2.8 x 10 6 Mp pc~ 3 , r = 0.22 pc and r/ = 1.75 
l|Schodel et al.l [20071 ). and obtain stellar disruption rates of 
~ 1.2 x lO _4 M0yr _1 . Our results presented in Table [2] give 
around 50 % higher rates in axisymmetric models, which 
lead to final BH masses in average 20 % higher with respect 
to spherically symmetric systems. This factor, would change 
the rates to ~ 1.8 x lO~ 4 A'/0yr -1 . Integrated over the age of 
the universe, the mass gain is of the order of the BH mass. 
It means, that relaxation processes might play an important 
role in the growth of the galactic center BH. Moreover, for a 
galaxy core of M ~ 10 9 Mq and a relaxation time of 10 Gyr, 
our models give us peak rates of the order of ~ lO -4 M0yr -1 . 
These rates are comparable to the accretion rates of some 
power law galaxies found by I Wang fc MerrittJ l|2004l . 

• Evolution of initially rotating systems with black holes 
affects substantially the orbital distribution of stars, spe- 
cially in the regions inside the BH influence radius. We have 
found that rotation in relaxation time scales can not be ne- 
glected, and it triggers higher disruption rates with an ex- 
cess of rotating stars. Central rotation has been detected 
in relaxation time scales in the zone of influence of the BH 
(fcrit < r < r n ), being at this time r n ~ r hm (0). The origi- 
nal frot.max moves inwards, in a region where the BH poten- 
tial dominates the stellar environment. For comparison, in 
systems without BH, thus without post-collapse evolution, 
dynamical insta bilities cause that u ro t,max moves outwards 
from the center (|Kim et al.l 12008! ). In the central profile of 
the parameter v Iot /a this maximum is lowered, mainly due 
to the cusp in the central velocity dispersion, but is not neg- 
ligible. 

The presence of different stellar populations is expected 
to enhance this effect in the central regions, as showed by 
iKim. Lee, fc Spurzeml (|2004T ). More massive stars segregate 
to the center in time scales shorter than a relaxation time, 
and they can rotate faster. Our currently investigations of 
multi-mass axisymmetric cores with stellar evolution and 
BHs, aim to obtain more detailed measurements of different 
stellar populations in the center, which could be detectable 
by observations. Another task, would be the treatment of a 
binary black hole (BBH), which can in a similar way, leads 
to a more efficient support in the development of rotation 



in its zone of influence (|Berczik et al.|[2~006l ; iBerentzen et al.l 
l2009h . 

More realistic NBody simulations by using a = 1 and 
higher particle numbers up to N ~ 10 6 or more, are never- 
theless necessary and still challenging to perform, specially 
in relaxation time scales, but will be possible in the near 
future. The advantage of using direct NBody models, to- 
gether with computationally faster Fokker-Planck realiza- 
tions makes possible to study the evolution of kinemati- 
cal and structural parameters in more detail, which can 
complement and test observational measurements. Observa- 
tional studies of 'collisional' galactic nuclei embedding mas- 
sive black holes, can be compared to evolutionary models to 
elucidate theoretical predictions and have a better under- 
standing of galaxy evolution. 
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